Differential equations are the bread and butter of physics. The purists would say that theoretical modelling (that is, arriving at the differential equations through physical reasoning) is the physicists' raison d'être. Perhaps they were right some decades ago. Today, the bulk of physics (and most of science at that) is application of existing models to ever more complicated systems. This could mean solving non linear partial differential equations in arbitrary geometries. It is clear that the subset of such problems that can be attacked with only pencil and paper is a very narrow one. The world is simply too complicated to be understood analytically.
Not all is lost though. We are the happy offsprings of the computer era. Many of the problems that shows us no mercy when armed with only pencil and paper offer little resistance to well coded numerical schemes. Our purpose here is to study some of the simplest and most robust of these. And then hopefully we shall also learn a thing or two about physics as well.
Let us consider first the simplest problem: that of a first order ODE in a n-dimensional space.
\begin{equation} \frac{d\mathbf{x}}{dt} = \mathbf{g}(\mathbf{x}, t) \end{equation}where $\mathbf{x}$ and $\mathbf{g}$ are written in boldface to highlight their vector nature:
\begin{equation} \mathbf{x}(t) = \big(x^0(t), x^1(t),...,x^{m-1}(t)\big), \qquad \mathbf{g}(\mathbf{x}(t), t) = \Big(g^0\big(\mathbf{x}(t),t\big), g^1\big(\mathbf{x}(t),t\big), ...,g^{m-1}\big(\mathbf{x}(t),t\big)\Big). \end{equation}Although it might seem at first sight that restricting ourselves to first order differential equations is too strong an impediment, the problem we are tackling here is actually very broad. Just consider how Hamilton's trick takes us from a system of second order differential equations in spacetime to a doubled system of first order differential equations in phase space. Such redefinition of higher derivatives into new dynamical variables is a useful tool that we shall illustrate soon.
At first, we shall look at this problem as an initial-value problem, meaning that our first order ODE is complemented by the initial condition
\begin{equation} \mathbf{x}(t_0) = \mathbf{x}_0. \end{equation}Assume we are given an array of instants of time $t$ and another one with the corresponding values of position $x$. We do not bother here to consider $x$ a vector, since the following discussion can be trivially extended to that case. We thus have
\begin{equation} x(t_0) = x_0, \quad x(t_1) = x_1, \quad..., \quad x(t_{l-1}) = x_{l-1} \end{equation}where the length of both arrays is $l$. Computing discrete approximations of the differential equation is easy peasy:
\begin{equation} \frac{dx}{dt}\big|_i = g_i \approx \frac{x_{i+1} - x_i}{t_{i+1} - t_i} \end{equation}where we have defined another array that stores the values of $g$
\begin{equation} g(x_0,t_0) = g_0,\quad g(x_1,t_1) = g_1,\quad ...,\quad g(x_{l-1}, t_{l-1}) = g_{l-1}. \end{equation}Since our goal is to numerically solve the ODE, we use the last equation to retrieve the values of $x_i$ at different times $t_i$. Thus we write
\begin{equation} x_{i+1} \approx x_i+ g_i(t_{i+1} - t_i). \end{equation}Let us implement the Euler method together with Hamilton's trick to solve the harmonic oscillator numerically. The harmonic oscillator ODE is
\begin{equation} \frac{d^2x}{dt^2} + \omega^2x = 0 \end{equation}where $\omega$ is the angular frequency. Here our dynamic variables are the position of the particle $x$ and its momentum $p = \frac{dx}{dt} = \dot{x}$. Thus the equation above can be written as a system of two ODES,
\begin{equation} \dot{x} = p\\ \dot{p} = -\omega^2 x. \end{equation}We can write the harmonic oscillator ODE as
\begin{equation} \frac{d}{dt}(x,p) = (p, -\omega^2 x). \end{equation}Now applying Euler's method, we get
\begin{equation} x_{i+1} \approx x_i + p_i(t_{i+1} - t_i)\\ p_{i+1} \approx p_i -\omega^2(t_{i+1} - t_i) \end{equation}The code written below is incomplete. Fill in the missing steps.
In [ ]:
After solving the exercise above, it should be clear that the Euler method isn't actually the best way to solve and EDO numerically. The reason is that it uses first order approximations to the derivatives. According to Taylor's theorem, we can approximate a $k+1$-differentiable function by it's $k$-th order Taylor polynomial,
\begin{equation} f(t) = f(t_0) + f'(t_0)(t-t_0) + \frac{f''(t_0)}{2!}(t-t_0)^2 + ... + \frac{f(t_0)^{(k)}}{k!}(t-t_0)^k + \frac{f^{(k+1)}(\xi)}{(k+1)!}(\xi-t_0)^{k+1} \end{equation}where $\xi$ lies between $x$ and $t_0$. We thus expect errors of order $O(\tau^{k+1})$, where $\tau\approx t-t_0$. Thus after $n$ steps we can expect an accumulated error of the order $nO(\tau^{k+1})$. Since $k = 1$ for the Euler method, we might quickly loose track of the real solution unless $\tau$ is very small. This is why the Euler method is actually almost never used.
Let us try to find other computational schemes based on the analytical solution,
\begin{equation} x(t) = x_0 + \int_{t_0}^tg(x,t')dt' \end{equation}The idea is clear: since we know $x_0$ and $g(x,t)$, in principle all we need to do is to use a quadrature scheme to compute the integral numericaly. Not so fast though. When we discretize this equation, we get
\begin{equation} x_{i+j} = x_i + \int_{t_i}^{t_{i+j}}g(x,t)dt \end{equation}and indeed the Euler method is trivially recovered by taking $j = 1$ and substituting $g(x,t) = g_i$. We thus see where the weakness of Euler's method stems from: it is based on the poorest of all quadrature schemes! But now it is also clear that we have a problem. Suppose we want to use the trapezoidal rule to compute the integral above. Although not very sofisticated, it is already an improvement: instead of taking $g(x,t) = g_i$, we now approximate $g$ as $\frac{g_{i+1} + g_i}{2}$, so we end up with
\begin{equation} x_{i+1} = x_i + \frac{g_i +g_{i+1}}{2}\tau. \end{equation}But $g_{i+1} = g(x_{i+1},t_{i+1})$, and we don't know $x_{i+1}$ yet. One way to bypass this problem is to use a fixed-point iteration, which goes by the name of Picard method among numerical analysis enthusiasts. It usually works pretty well, and consists in performing the recursion
\begin{equation} x_{i+1,0} = x_i\\ x_{i+1,j+1} = x_{i,j} + \frac{g_i + g_{i + 1,j}}{2}\tau. \end{equation}a certain number of times.
The logistic function, also known as the sigmoid curve, is defined as
\begin{equation} x(t) = \frac{L}{1+\text{e}^{-\omega(t-t_0)}}. \end{equation}It is ubiquitous in applied mathematics, ranging from artificial neural networks to sociology. In physics, it is known as the Fermi-Dirac statistics.
One intuitive way to view this function is to derive it as the solution of the first order non-linear ODE
\begin{equation} \frac{dx}{dt} = x(1-x) \end{equation}with boundary condition $x(0) = x_0 = \frac{1}{2}$. The above ODE can be seen as a model for population growth in which the relative change in population is proportional to both the current population $x$ and the available resources $1-x$ in the environment.
The code below uses the Picard method to numerically solve the logistic differential equation. Fill in the remaining steps.
In [ ]:
We now introduce the industry standard Runge-Kutta method. Its advantage is that it uses a clever strategy to improve accuracy through a series of first order approximations of the function $x(t)$. As we shall see, this is accomplished by evaluating $g(x,t)$ at several points, in opposition to evaluating higher order derivatives of $g$ at the same point. One of the nice things about the Runge-Kutta method is that it doesn't demand extra data to work. This is a big plus, since in most dynamical problems we are only given a single initial condition.
Again our starting point will be Taylor's expansion. For the sake of clarity and brevity, let us use the Runge-Kutta method of order 2 as a proxy. It starts from the second order Taylor expansion,
\begin{equation} x(t + \tau) = x(t) + \tau x'(t) + \frac{1}{2}\tau^2x''(t) + O(\tau^3). \end{equation}Now remember that we are given the following data, that actually specify the initial value problem we want to solve:
\begin{equation} x'(t) = g(x,t), \qquad x(t_0) = x_0. \end{equation}The idea then is to use a Taylor expansion of $g(x,t)$ to substitute $x'(t)$ and $x''(t)$. This leads us to the followinyg equation.
\begin{equation} x(t + \tau) = x + \tau g + \frac{1}{2}\tau^2(g_xg + g_t) + O(\tau^3), \end{equation}where $g_x(x,t)$ and $g_t(x,t)$ are the partial derivatives of $g$ with respect to $x$ and $t$, respectively.
Now observe the following. Instead of using the usual Taylor expansion of $x$ at a given order, we can try to make successive first order approximations. For example, the first order approximation for $x(t + \tau)$ is simply
\begin{equation} x(t + \tau) \approx x(t) + \alpha_1h_1 \end{equation}with $\alpha_1 = 1$ and $h_1 = \tau g = \tau \frac{dx}{dt}$. Note that we are using the quantity $\tau g$ a "standard" first order variation in $x$. If we stopped here, we would just recover Euler's method. So let us carry on. To the next order we write
\begin{equation} x(t + \tau) \approx x(t) + \alpha_1h_1 + \alpha_2h_2 \end{equation}where
\begin{equation} h_1 = \tau g(x,t) \qquad\text{ and} \end{equation}\begin{equation} h_2 = \tau g(x + \nu_{21}h_1, t + \nu_{21}\tau). \end{equation}If we now Taylor expand $h_2$ we end up with
\begin{equation} x(t+\tau) \approx x(t) + (\alpha_1 + \alpha_2)g\tau + \alpha_2\nu_{21}(g_xg+g_t)\tau^2 \end{equation}which must be compared with our original Taylor Expansion in equation (???). Such comparison leads us to the following equations
\begin{equation} \alpha_1 + \alpha_2 = 1 \end{equation}\begin{equation} \alpha_2\nu_{21} = \frac{1}{2} \end{equation}An interesting feature of the Runge-Kutta method is that the choice of parameters is not unique, since we have more parameters than constraints. A particularly symmetric solution is $\nu_{21} = 1$, $\alpha_1 = \alpha_2 = \frac{1}{2}$, for which $x(t+\tau)$ is approximated as
\begin{equation} x(t + \tau) \approx x(t) + \frac{\tau}{2}\Big(g\big(x(t),t\big) + g\big(x(t)+\tau g(x,t), t + \tau\big)\Big). \end{equation}As we can see, this choice corresponds to taking the approximated average of $x'(t)$ at $t$ and $t + \tau$ and using it in the first order Taylor formula. Also note that we have accomplished to substitute the evaluation of partial derivatives of $g$ at the same point by the evaluation of $g$ itself at two different points. This hallmark of the Runge-Kutta method makes its implementation very easy. To illustrate this point, consider the widely used Runge-Kutta method of order 4, which is given by the following set of equations.
\begin{equation} x(t + \tau) = x(t) + \frac{1}{6}(c_1 + 2c_2 + 2c_3 + c_4) \end{equation}where
\begin{equation} c_1 = \tau g(x,t), \end{equation}\begin{equation} c_2 = \tau g\big(x + \frac{c_1}{2}, t + \frac{\tau}{2}\big), \end{equation}\begin{equation} c_3 = \tau g\big(x + \frac{c_2}{2}, t + \frac{\tau}{2}\big) \end{equation}\begin{equation} c_4 = \tau g(x + c_3, t + \tau). \end{equation}You can check that this choice of parameters leads to a correct Taylor expansion of $x(t)$ to the fourth order.
In [ ]: